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ABSTRACT 

We develop a Monte-Carlo technique based on L.B. Lucy's indivisible photon pack- 
ets method to calculate X-ray continuum spectra of comptonized thermal plasma in 
arbitrary geometry and apply it to describe the broadband X-ray continuum of the 
galactic superaccreting microquasar SS433 observed by INTEGRAL. A physical model 
of the X-ray emitting region is proposed that includes thermal emission from the ac- 
cretion disk, jets and hot corona where the photons of different origin are comptonized. 
From comparison with INTEGRAL observations, we estimate physical parameters of 
the complex X-ray emitting region in SS433 and present model spectra for different 
viewing angles of the object. 
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1 INTRODUCTION. 

SS433 is a peculiar m assive X-ray b inary system with pre- 
cessing relativistic jets ('Margon"l984') . The system is located 
at a distance of about 5.2 kpc ( Lockman et al.ll2007h nearly 
in the galactic plane. The optical companion V1343 Aql was 
first id entified in the survey of st a rs exh ibiting Ha emis- 
sion of IStephenson and SanduleakI (|l977h . It is one of the 
brightest stars in the Galaxy, the bolometric luminosity of 
the object assuming isotrop ic radiation is Lboi ~ 10^" erg/s 
l|Cherepashchuk et al.|[l98a ). SS433 i s a close binary system 
with an orbital period of 13.1 days l|CherepashchukHl98ll) . 
The uniqueness of this source comes from the presence of 
narrow oppositely directed subrelativistic jets emitting red 
and blue-shifted periodically variable lines. The commonly 
accepted model of SS433 suggests a continuous regime of 
supercritical accretion of gas onto the relativistic star. The 
optical star fills its critical Roche lobe supplying powerful 
and almost continuous flow of gas onto th e relativistic star 
at a rate of ~ W^MQ/yr (|Fabrikall2004l 'l. A supercritical 
accretion disk forms together with narrow jets of gas prop- 
agating perpendicular to the disk plane from the central 
regions of the disk and having a relativistic speed of 0.26c. 

Presently, SS433 is recognized as a galactic microquasar 
with precessing supercritical accretion disk around a rela- 
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tivistic compact object, and has been extensively investi- 
gated in the optical, radio and X-ray ranges (fo r a compre- 
hensive review and references see lFabrikal|2004l )). Masses of 
both components were obtai ned from the analysis of the op- 
tical light curves of SS433 (|Antokhina and CherepashchukI 
1987) and of the hard X-ray eclipses observed by INTE- 
GRAL (Cherepashchuk et al. 2008), so the relativistic object 
in SS433 is most certainly a black hole. Despite a significant 
progress in understanding the nature of SS433, the mecha- 
nism of collimation and acceleration of matter in jets to the 
relativistic velocity is still unclear. 

In the present paper we shall concentrate on the broad- 
band X-ray spectrum of SS433. In th e last years t he 



source was intensively o b serve d by RXTE |Revnivtsev et al 
'2004 iFiUppova et al.l bOOel) . XMM ferinkmann et al 



1200'n). Chandra (iMarshall et al I |2005|) and INTEGRAL 



I Cherepashchuk et al.ll2003l , |2004 12OO5I |2007| ) . INTEGRAL 
observations of SS433 discovered a hard X-ray spectrum 
from the source up to 100 keV visible over the whole preces- 
sional period ~ 162 days. This suggests the presence of an 
extended hot r egion in the central parts of the supercr itical 
accretion disk (jCherepashchuk et al.|[2005l . l2007l , l2008h . 

The initial analysis of the broadban d X-ray contin- 
uum observed by INTEGRAL (jCherepashchuk et al.ll2005l . 

allowed for the first time to estimate the tempera- 
ture and Thomson optical depth of the assumed hot scatter- 
ing corona above the accretion disk in SS433. However, the 
spectral analysis made use of simplifying assumptions about 
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Figure 1. Schematic picture of the model adopted for SS433 

the corona geometry (mainly because of utilizing the stan- 
dard models incorporated in the XSPEC tool). The ther- 
mal emission from the X-ray jet was treated independently 
and it was possible to estimate the mass loss rate in the 
jet. The joint modeling of hard X-ray eclipses and preces- 
sional variability allowed the estimation of the size of hot 
corona in the framework of the adopted g eometrical model 
of the source (ICherepashc huk et al.l l2007. 200^ It is well 
established from observations (see iFabrika (2004 ) for more 
detail) and es pecially from recent studies o f X-ray variabil- 
ity of SS433 (|Revnivtsev et al.ll2004l . |2006| ) that X-ray and 
optical variabilities are correlated on time scales less than 
10^ s. The power spectrum analysis suggests that there is an 
extended accretion disk with size similar to that in Cyg X- 
1 in the soft state l|Titarchuk et al.ll2007l l. The strong slow 
wind outflow from the supercritical accretion disk forms a 
photosphere, and there is a funnel filled with hot plasma 
around the jet. Clearly, to obtain the correct physical pic- 
ture, more accurate modeling of the emission processes from 
the jet, accretion disk and hot corona is needed. 

The purpose of this paper is to model the observed 
broad-band X-ray continuum of SS433 using accurate treat- 
ment of emission processes in the jet and the scattering hot 
corona above the accretion disk. To calculate the broad- 
band X-ray spectrum, we shall use the Monte-Carlo method 
(see the Appendix for a more detailed description). The first 
reason for our choice is connected with the geometry of the 
system containing the accretion disk, the corona and the 
jet. This is a rather complicated geometry, it has azimuthal 
symmetry, but for the outgoing radiation nevertheless it is 
impossible to obtain any analytical solutions, especially with 



account of the temperature decrease along the jet. The sec- 
ond reason is connected with physical processes, in partic- 
ular the Compton scattering of thermal photons from disk 
and jet by free electrons in the hot corona. The accurate 
analytical treatment of this process in the radiation trans- 
fer equation makes it impossible to solve it. For the same 
reasons, the computational algorithm for numerical solu- 
tion of the transfer equation will be very sophisticated and 
will demand much computational time. An advantage of the 
Monte-Carlo method is that the computational algorithm 
is rather simple. It consists roughly of one procedure re- 
peated many times, and quite a small amount of compu- 
tational time is required to calculate the emergent spectra 
with good accuracy. Moreover, when using finite-difference 
methods, computational errors tend to increase with time 
as they are accumulated at certain steps of the algorithm. 
In the Monte-Carlo case the picture is completely different: 
the precision of the result increases with computational time 
because of larger number of individual tests and thus better 
statistics, so waiting longer gives better result. 



2 GEOMETRICAL MODEL OF SS433 

The schematic picture of the source is presented in Fig[l] 
The coUimated, o pposi tely directed j e ts we re proposed by 
iFabian and Reed (|l979l ') and iMilgromI l| 19791 ) to explain the 
observed moving emission lines at unusual wavelengths. 
These lines were interpreted as Doppler-shifted Balmer and 
He I lines (|Margon et al.lll979h . The morphology of SS433 
and the surrounding nebula from radio observations and the 
analysis of its temporal behavior provides strong support to 
the idea of ballistically coasting matter ejected from the twin 
jets, which jus tifies the so-called "kine matic model" of the 
jet precession (|Abell and Marg"onlll979l ) proposed to explain 
the periodic "movement " of the opti c al lin es. The X-ray 
imaging observations of ISeward et all l|l980l) also support 
the jet model. The jets must be supplied with matter from 
some source, so a hypothesis of a compact star in a close 
binary system was proposed. The discovery of the 13-day 
periodicity in the ra dial velocity of the "st ationary" SS433 
emission-line system (jCrampton et al.lll980t l proved this pic- 
ture. The donor star fills its critical Roche lobe providing 
powerful flow onto the compact object. Thus, an accretion 
disk forms with a pair of narrow jets propagating perpen- 
dicular to the disk plane. The presence of the accretion disk 
follows also from the interpretation of the light curve that 
involves mutual eclipses of the massive carly-typc normal 
star and the extended accretion disk ((Chcrcpashch uk 198l|; 
lAntokhina and Cherepashchuklll987l ). As was already men- 
tioned in the Introduction, the presence of a hot corona 
above the inner regions of the accretion disk was sug- 
gested to explain the hard X-ray spectrum of SS433 ob- 
served by INTEGRAL over the entire precession period 
ICherepashchuk et akllioosl . |2007| . l2008h . So the considered 
model of SS433 including the accretion disk with hot corona 
above its central regions and two jets perpendicular to the 
disk plane is well justifled by observations. 
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2.1 Supercritical accretion disk 

In SS433, the optical star fills the Roche lobe and supplies 
matter to the accretion disk at a rate of ~ 10~*MQ/yr. 
The accretion disk radius is assumed to be limited by the 
size of the Roche lobe and thus for Mi,h = 7Mp) and q = 
Mbh/Mopt = 0.3 (jAntokhina and Chereoashchukl 1 1987^ 1 is 
~ 10^'^ cm. The most emission from the disk is generated in 
its central regions, so the exact value of the accretion disk ra- 
dius is of small importance. Far away from the central parts, 
the s tandard accretion disk theory ((Shakura and SvunvaevI 
1 19731 ) can be used to estimate its geometrical thickness. For 
the disk height h growing linearly with radius the disk open- 
ing angle Odisk can be directly found from the relation: 



h = r tan ■ 



(1) 



We shall consider an optically thick disk, so in our calcula- 
tions this conical surface will be used as an opaque screen for 
incident photons. For reasonable input parameters discussed 
below, the disk opening angle is 9disk = 4° . 



2.2 Thermal X-ray jets 

The evidence for X-ray emitting jets in SS433 was found 
in X-ray observations of the source. The EXOSAT satel- 
lite detected Doppler shifted iron lines moving periodically 
in accordance with the precessional Doppler curve of the 
source (Watson et al,.1986. ), suggesting the X-ray emission 
origin from the base of jets. Later observ ations by Ginga 
l|Kawai et al.ll 19891: iBrinkmann et aLlliggil ) and deep ASCA 
observations i Kotani et al.l T994I . 19961 ) discovered many 
pairs of Doppler-shifted emission lines originating in the blue 
and r ed jets moving with a radial velocity of 0.26c (|Margonl 
1 19841 ). The soft X-ray emission is fitted by a combination 
of a thermal bremsstrahlung continuum and emission lines 
from a multi-temperature optically thin plasma in the coni- 
cal jets with length Lj ^ 10^^ cm and a total X-ray luminos- 



ity of ^ 3 X 10^^ erg/s llKotanilll996 ). Mo re recent Chandra 
jMarshall et al.ll2005l: iNarniki et al.lf2003l ') and XMM obser- 
vations ( Brinkmann et al.l 20051 ) confirmed the basic prop- 
erties of X-ray jets and put constraints on their geometrical 
and physical parameters. The much higher kinetic luminos- 
ity of jets ~ 10'^'"' — 10*" erg/s is derived from thei r inter action 
with the surrounding nebula W50 (see iFabrikal (|2004l ') for a 
detailed discussion and further references) . We shall assume 
X-ray jets to have a conical shape, as inferred from the coin- 
cidence of the jet opening angle 6j in the optical and X-rays 
l|Fabrikall2003 ). with the jet opening angle Oj = 1.2°. 



2.3 Hot corona 

The presence of a hot corona above the central parts 
of the accretion disk in SS433 has been suggested by 
hard X-ray observations of the sourc e by INTEGRAL 
ICherepashchuk et alllioosl . |2007| , |2008| ). The main indica- 
tion is the unchanged from of the hard X-ray spectrum at 
different precession phases. If only thermal emission from 
cooling jets were responsible for the observed X-rays, one 
would expect a significant softening of the spectrum at the 
precessional phase r3 (disk edge-on), where the innermost 
hottest parts of X-ray jets are screened by the accretion 
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Figure 2. The computational domain. 
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■ edge. The observed primary X-ray eclipse light curve 
and precessional hard X-ray variability cannot be success- 
fully described by thermal jets without a broad hot structure 
above the central parts of the di sk which is partially eclipsed 
at th e r3 precessional phase (jCherepashchuk et al.l l2005l . 
|2008| ). The hard X-ray spectrum at the maximum opening 
disk phase can be modeled as scattered thermal radiation 
from jets by an isothermal hot corona {Tcor ~ 18 — 30 keV) 
with spherical or slab geometry using different standard 
spe ctral tools from the XSCPE C package (see, for exam- 
ple, [Chereeashc^kelaD (|2007l ) for the ComPS model fit). 
The corona around jets can be heated up to such temper- 
atures, for example, by the interaction of the precessing 
jets with wind outfiow from the supercritical accretion disk 
iBegel man et al] (|2006l ). Motivated by additional geometri- 
cal constraints obtained from the analysis of the joint pre- 
cessional and eclipsing X-ray variability of SS433 which sug- 
gests a fairly broad opening angle of the jet funnel, we shall 
assume a spherical shape of the corona. The parameters of 
the model are shown in Fig[T] 



3 PHYSICAL MODEL OF THE X-RAY 
EMITTING REGIONS IN SS433 

We consider the accretion disk with a spherical corona above 
its inner regions, and two jets propagating perpendicular 
to the disk plane. This picture is symmetric with respect 
to the disk plane, so we shall consider only its upper half 
(Fig[l|. The system is axially symmetric with respect to the 
jet axis, so radial profiles of physical quantities fully charac- 
terize this 3-dimensional structure. The computational do- 
main (see Fig[2]) is bounded by the surface of the accretion 
disk (from bottom) and a sphere with radius equal to the jet 
height. The inner accretion disk is assumed to have constant 
hjr ratio, so the disk surface is inclined by 2° to the equato- 
rial plane. The conical jet, disk and corona begin at the coor- 
dinate origin. Both jet and corona are assumed to consist of 
fully ionized hydrogen plasm a. Guided by observations of the 
accre tion wind photosphere (|Fabrikall2004l : iRevnivtsev et al.l 
l2004h . we set the visible jet base at ro — 10^^ cm and as- 
sume the inner spherical boundary of the corona to lie at 
the same radius, as shown in Fig[2l 
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3.1 The jet 

We fix the jet opening angle at 1.2 degrees. The temperature 
at the base of the jet is assumed to be the same as that of 
the corona and is set to be To = IQkeV , the best-fit value 
found in our calculations. 



T = To 



4/3 



2/3 



(2) 



In our simulations we shall assume that the jet tempera- 
ture changes with height according to the law correspond- 
ing to the jet cooling due to adiabati c expansion of ideal 
monoa tomic gas ([2)| . In the paper of iKoval' and Shakural 
l| 19891 ) an analytic solution for the jet temperature profile 
was obtained taking into account both adiabatic and ap- 
proximate treatment of free-free emission cooling. For the 
constant jet opening angle the solution obtained by Koval' 
& Shakura is determined by the mass rate in the jet, Mj. If 

rapid cooling of the jet plasma occurs. Here TZ is the uni- 
versal gas constant and To — T^or is the temperature at 
the base of the jet. For jet parameters considered by us 
Miim — 5.3 • 10^* g/s. From observations, the ki netic lu- 
minos ity of the jet is about Lkin ~ iq-^^-i-o gj-g^g jFabrikal 
so the lower limit for the mass loss rate in the jet 
is Mj = 3 • 10^^ g/s, which is six times greater than Mum- 
For Mj/Miim < 1 the jet cooling occurs almost as in the 
adiabatic case, but in the opposite case, M /Mum > 1, the 
jet begins to cool rapidly with radius and its total X-ray 
luminosity should be so small that it would be impossible 
to observe it. The jet can have a larger opening angle so 
that the value of Mum increases and may exceed Mj , which 
will prevent the jet from rapid cooling. There could be also 
some processes (Compton heating by hard quanta from the 
corona, interaction with the wind from the disk during the 
jet precession, etc.) that additionally heat the jet and pre- 
vent it from rapid cooling. The accurate treatment of the jet 
heating has not been done yet, so we assume the adiabatic 
temperature distribution along the jet. 

The density dependence along the jet follows from 
the continuity equation. The jet mass outfiow rate Mj = 
nnipVS , where n is the number density of protons at r, rrip 
is the proton mass, V = 0.26c and S = 7rr^ tan 6jet/2 is the 
jet cross-section. So we can write 



n = no 
where 



M. 



no 



jet 



TvnipVrQ 



tan 



Ujet 



(3) 



(4) 



parts of the disk. When the disk luminosity approaches 
the value of ~ 0.6 Eddington limit, the radiation force 
exceeds the gravitational one and the outfiow from the 
central parts of the disk sets in. The convection motions 
in the disk generate acoustic fiux, which dissipates in the 
upper layers of the disk and heats them. As a result, a hot 
(T ~ 10^ K) rarefied corona is formed above the central 
parts of the accretion disk. The second me chanism with 
accou nt for magnetic field was studied by iGaleev et al.l 
It is connected with magnetic field amplification 
in the disk. The magnetic field generated by differential 
rotation of conductive media is amplified due to convective 
motions in the direction perpendicular to the disk plane. 
The amplified field can attain strength comparable to 
equipartition value and magnetic flux will emerge from the 
disk, leading to formation of accretion disk corona, which 
consists of many magnetic loops, the flelds of which provide 
an energy source for the corona heating. The corona is hot 
and has a low density compared to the accretion disk. In 
the present paper we assume a spherical isothermal corona 
with temperature 19keV = 2.2 • 10* K. For simplicity, the 
radial density dependence in the corona is assumed to obey 
the same law as in the jet but with the different initial 
value at r = lO^'^ cm. Then the optical depth for Thomson 
scattering in the corona is 



Tcor = (TTno 



(5) 



The radius of the corona can be expressed in terms of its 
optical depth as 



fcor 



ro/(l 



aTncoi-ro 



(6) 



near = (4.3 - 4.8) ■ 10 cm" 
<tt ~ 6.65 • 10~^^cm 



The analysis of the observed wide primary orbital 
hard X-ray eclipse and the precessional variability 
I Cherepashchuk et al.|[2007l . l2008h suggests large visible size 
of the hot corona comparable to that of the accretion disk, 
but it is not absolutely certain. We have fitted the observed 
X-ray spectra with two models of corona: 1) the corona 
is not eclipsed by the precessing accretion disk and has 
Tcor ~ 6.4 ■ 10^^ cm, Tcor ~ 0.24; 2) the corona is fully 
eclipsed with rcor ~ 2.7 • 10^^ cm, Tcor = 0.20. Both models 
can reproduce the observed X-ray spectrum, but the first 
model gives a slightly better fit. We shall discuss in more 
detail the obtained results below. 



3.2 The corona 

The existence of hot coronas above highly lu- 
minous accretion disks w as fir st stu died by 
iBisnovatvi-Kogan and Bhnnikovl l| 19761 . Il977l ). where 
two mechanisms responsible for corona formation have 
been considered. The first one is connected with disk 
radiation pressure. The particle moving in the vicinity 
of a radiating disk suffers the action of radiation from 
the whole disk. Thus, the radiation force depends on the 
position of the particle and is maximal at the central 



3.3 The accretion disk 

Basic parameters of the standard accretion disk model in- 
clude the mass of the compact object M, mass accretion 
rate through the disk M, and the angular momentum trans- 
port efficiency parameter a. We are interested in the 'effec- 
tive' mass accretion rate through the disk and it is different 
from the mass rate supplying by the optical star through 
the inner Lagrangian point. The jet kinetic power is ap- 
proximately equal to the Eddington luminosity for the cen- 
tral black hole mass Mj,^ = 7Mq. This power corresponds 
to the critical accretion mass rate through the disk Mcr- 
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Figure 3. X-ray continuum of SS433 from 3 to 90 keV 



As the same mass rate should reach the central source 
to be ultimately consumed by the black hole, the 'effec- 
tive' mass flow rate through the disk is M ~ 2Mcr- Tak- 
ing M = 7Mq, M « 2Mcr and a = 0.1 we obtain that 
in the outer gas-dominated region with free-free opacity 
l|Shakura and SvunvaevlligTST ). which lies beyond the radius 



rbc = 4-10 cm for the adopted parameters, 



h = 6.1- lO^a 



-1/10 



R 



-1/2 



M 
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9/8 



(7) 



At the accretion disk radius Rdisk = 1.5 ■ 10^^ cm the half- 
thickness of the disk is 5.5-10^" cm, so the angle between the 
disk surface and the equatorial plane in ([T| is ddinhj'i, ~ 2° . 
To calculate the emergent photon spectrum from the corona, 
we shall assume that all accretion disk radiation is emitted 
only in the central region, from a sphere with radius ro, 
while the outer regions of the accretion disk play the role of 
an opaque screen for photons (Fig[2]). Accretion disk photons 
are distributed uniformly across the sphere with radius ro 
and initially move strictly in radial direction. 



4 RESULTS OF SIMULATIONS 
4.1 Observational data 

For spectral analysis, we have made use INTEGRAL obser- 
vations of SS433 in May 2003 obtained simultaneously by 
IBIS/ISGRI (20-100 keV) and JE M-X (3-20 keV) telescopes 
l|Cherepashchuk et al.|[2003l . |2005| ) . The source was observed 
near the precessional phase (the T3 moment) when the ac- 
cretion disk is maximally open and the flux from the source 
is also maximal. The angle between the line of sight and 
the blue jet (i.e. the jet pointing to the observer) at this 
phase is ~ 60'^. The Doppler redshift at this precessional 
phase is z ~ —0.1. Spectral lines near 7 keV were ignored. 
The X-ray data were processed using the standard OSA 5.1 



Figure 4. The comparison of the Monte-Carlo simulated X-ray 
continuum of SS433 with observations, dabs = 60", only the jet 
and corona emission are considered. The radius of spherically- 
symmetric corona is Tear = 6.4 • 10^^ cm, the mass loss rate in 
the jet is Mj^t = 4 • 10-*^" g/s. The viewing angle of the system is 
dobs = 60". 




Figure 5. The same as in Fig. |4] for smaller corona with rear = 
2.7 ■ lO'^'^ cm, and Mj^t = 3 • 10^^ g/s, viewed by the same angle 
0obs = 600. 



INTEGRAL softwar e developed by the INTEGRAL science 
data center (ISDG, http://isdc.unige.ch| l|Courvoisier et all 
l2003l j). The resulting 3-90 keV X-ray spectrum is shown in 
FigEl 

A model fit of the observed spectrum is accepted by 
achieving the minimal reduced x iXred = X /N),N_=_n — 
rijtar i s the number of the degrees of freedom, see IZombec^ 
(|2006t ). 



4.2 Bremsstrahlung comptonization model 

In this section we present the results of X-ray spectral mod- 
eling neglecting the radiation from the disk. The distance to 
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Parameter 


Eclipsing corona 


Non-eclipsing corona 


'"con -10^^ cm 


2.7 


6.4 




0.24 


0.20 




4.8 


4.3 


Mj,-10l9 g/s 


3 


4 


Lkin, -10^^ erg/s 


0.9 


1.2 




21.40 


23.34 


CL 


91.6% 


95.4% 



Table 1. Physical parameters of two corona models 



SS433 is fixed at d = 5 Itpc = 1.5-10^'' cm. The emerging X- 
ray spectrum in the 3-90 IteV range is assumed to be formed 
due to free-free emission of corona and jet with account for 
Compton scattering. After many trials, the following best-fit 
parameters of the model have been obtained (quantitative 
characteristics of the fitting are discussed in the next para- 
graph). The temperature of the corona and the jet base is 
Tcor ~ 19 keV. Fig|4] shows the comparison of the computa- 
tional results with the observed broadband X-ray spectrum 
of SS433 for a corona with rcor ~ 6.4 ■ 10^^ cm, Tcor ~ 0.24, 
ricor = 4.3- lO^'^ cm~'^ which is not fully eclipsed by the disk. 
The mass loss rate in the jet is Mj — 4 • 10^^ g/s, which cor- 
responds to the kinetic luminosity of Lkin ~ 1-2 ■ 10"^^ erg/s. 
The next Fig[5] shows the results of simulations for a fully 
eclipsing corona with radius rcor = 2.7 • 10^^ cm, the 
Thomson optical depth Tcor = 0-2, and the number den- 
sity ricor = 4.8 • lO'^^cm"''. The mass loss rate in the jet is 
Mj = 3- 10^^ g/s corresponding to the jet kinetic luminosity 
Lk^n = 9 • 10^^ erg/s. 

It is clearly seen that in both cases the model can de- 
scribe the observed X-ray continuum. The data includes 39 
experimental points. There are 5 fitting parameters: the jet 
mass loss rate Mj, the coronal temperature Tcor, number 
density at ro, Ucorfl, the Thomson optical depth Tcor, and 
the initial radius ro, so the number of the degrees of free- 
dom is Ndof = 34. For the model with non-eclipsing corona 
= 21.4 shown in Fig. g] the reduced xlcd = 0.630, and 
the confidence level of the fit at Ndof = 34 is « 95.4%. 
For the second model with fully eclipsing corona showmn 
in Fig. [S] = 23.34 for 34 degrees of freedom, so the re- 
duced Xred ~ 0.687 and the corresponding confidence level 
of the fit is « 91.6%. Thus the model with non-eclipsing 
corona (FigU)) with smaller seems to be more realistic. 
All model parameters are listed in Table |4]2] 

The angular dependence of the emergent X-ray spec- 
trum is shown in Fig[S] The spectrum calculated for the 
viewing angle 82° is similar to 60°, the minimum viewing 
angle of the jet during its precessional motion. The weak 
dependence on the viewing angle may be due to the corona 
(which radiates almost isotropically) mainly contributing to 
the total emission at all angles. This figure also shows that 
the fiux from the system is at maximum for the viewing 
angle 90° when both jets and both parts of the corona con- 
tribute to the emergent spectrum. The observations, how- 
ever, suggest that the fiux from the source at this preces- 
sional phase (the moment Tl, disk edge-on) is minimal. This 
may indicate that the supercritical accretion disk should be 
much thicker geometrically th a n pre dicted by the standard 
model [Shakura and SvunvaevI l| 19731 ). Significant screening 
of the corona by a thick disk is needed to explain the 
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Figure 6. The two-component model spectra of SS433 calculated 
for different viewing angles: 60 degrees (the bold solid line), 82 
degrees (the dash-dotted line), 90 degrees (the dotted line). The 
emergent spectrum for 90 degrees was obtained by doubling the 
computational data because of equal contribution of both jets and 
coronas, rcor = 6.4 ■ 10^^ cm, Mjct = 4 • 10^^ g/s. 



observed precessional variability of SS433 in hard X-rays 
(| Cherepashchuk et alllioOTl . l2008l l . 

The calculated spectrum of the source does not vary 
significantly with the viewing angle in the allowed range 
dobs ~ 60° — 90°, but if we were looking straight along 
the jet {9obs = 0°), we whould ob s erve a com pletely differ- 
ent picture CFigrfll . iFabrika et al.l (|2004 l2006h argued that 
ultraluminous X-ray sources (ULX) are superaccretors like 
SS433 seen almost along the jet. This hypothesis is based on 
the following assumptions: the most central parts of the su- 
percritical accretion disk are surrounded by outflowing gas 
which is semi-opaque, so the 'edge-on' luminosity of such an 
object should be much smaller than the 'face-on' luminos- 
ity, when the hottest central regions are exposed directly. 
We do not consider the outflowing gas, which is not visible 
in observational data, so the angular dependence of X-ray 
luminosity is completely different, as demonstrated above. 
Moreover, if we were looking straight along the jet, the lu- 
minosity would be actually minimal, since along the axis 
the jet itself is optically thick for its own thermal emission. 
The thermal energy power emitted per unit volume in the 
jet decreases with distance as both density and temperature 
decrease, so the most jet power is released near the jet base. 
Thus, in order to see the hottest central regions of the ac- 
cretion disk, the system should be observed not along the 
jet axis, but slightly from aside, by viewing angles exceeding 
the jet opening angle. 



4.3 Influence of the accretion disk emission 

The two-component model described in the previous section 
gives satisfactory results, but the accretion disk emission 
could also affect the emergent X-ray continuum. To test this 
in our calculations, we included the emission from the inner 
parts of the accretion disk with the following spectrum. 
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Figure 7. The model speetrum of SS433 observed along the jet 
axis (the solid line) compared with spectrum at Bobs = SO'^ (the 
bold solid line). 



\ = ALdisk/lfc exp (l - ^) , > i^c 

""^^^ 

(Here Ldisk is the total luminosity of the accretion disk.) 

The value r^, which is interpreted as the inner ra- 
dius of accretion disk, was taken to be rsph/0.3, where 
Tsvh = GMM / Lcr is th e spherization radius (see 
IShakura and SvunvaevI (| 19731 )) at which the disk lumi- 
nosity equals the critical Eddington luminosity Lcr = 
A-KcGMih/Kes = 8.8 ■ 10^*^ er g/s. It was shown by 
iBisnovatvi-Kogan and Blinnikovl l| 19771 ) that the outflow 
from a supercritical accretion disk begins at Ld = Q.GLcr- 
Here we take the lower luminosity of the disk Ld = O.lLci-. 
This may be justified by analysis of complicated physical 
processes in the supercritical accretion disk, but the main 
reason for our choice is based on our simulations: at higher 
Ld an excess of soft X-rays appears, which is absent in ob- 
servations. For supercritical accretion the outflow can start 
from radius where the luminosity is still subcritical, since 
the disk should be additionally heated by radiation from 
hot corona. The best fit with the observed spectrum oc- 
curs when the luminosity is Q.lLcr- Using the standard disk 
accretion equations and Paczynski-Witta gravitational po- 
tential tfipw ~ GM/{r — rg), the effective temperature can 
be found to be 

r.Hr) = 2.05 . 10 (j^^^j [ — ) (8) 

Here nn = S rp = 2GM/(? and (/'r,ri„ is 

(|Bisnovatvi-Kogan|[l989l) 




Figure 8. The three-component model of SS443 spectrum. The 
model spectrum is almost identical to the two-component model 
shown in FigO Taor = 2.7 ■ 10" cm, Afjet = 3 ■ 10^^ g/s. 

For the assumed disk parameters we obtain from these for- 
mulas Te/ 6.0 • 10^5 K « 0.05 keV. 

The accretion disk emission is assumed to be emerging 
from the sphere with radius ro = lO'^^ cm with all photons 
initially moving radially. Inside this radius the heat pro- 
duction in the disk exceeds the local Eddington value driving 
the mass outflow. We assume that rc is approximately equal 
to the inner radius of the jet. We do not take into account 
the radiation from the region inside Tc, assuming that it is 
transformed into the jet energy and heating of the corona. 
The results of the three-component model simulations are 
presented in Fig|8]and[9] 

It is clear that generally the disk radiation can strongly 
change the emerging spectrum, but in the energy range 1- 
100 keV the spectral difference is insignificant because of 
a comparatively low effective temperature T^f ~ 0.05 keV 
of the disk emission. In principle, the spectrum at lower 
energies can be used to estimate the accretion disk physical 
parameters. 



5 DISCUSSION 

We have developed a Monte-Carlo code for modeling X-ray 
(3-90 keV) continuum of SS433 based on the realistic physi- 
cal model of the source, including the accretion disk, thermal 
jet and hot scattering corona. We found the model param- 
eters that provide good qualitative and quantitative agree- 
ment with INTEGRAL observations. The fitting procedure 
was as follows. 

The similarity of hard X-ray spectra of SS433 at dif- 
ferent phases of the precessional period, as well as the joint 
modeling of the orbital primary X-ray eclipse and preces- 
sional variability, suggest the presenc e of a hot corona that 
is not fully eclipsed by the donor star l|Cherepashchuk et al.l 
|2008| ) . This fact allows us to put the lower limit on the outer 
radius of the corona. 
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Figure 9. Comparison of two- and tliree-component emission 
models. The solid curve corresponds to the model with disk ra- 
diation, the dashed one — to the model without disk radiation. 
r^or = 6.4 ■ 10" cm, Mj^t = 4 ■ 10^^ g/s. 
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Figure 10. The results of the simulations for different corona 
temperatures: 17 keV (the dash-dotted line), 19 keV (the solid 
line), 21 keV (the dotted line). The best fit is found for T^or = 
19 keV, rcor = 0.24, rcor = 6.4 ■ lO" cm, Mj et = 4 ■ lO^'' g/s. 

First, we adjust the shape of the spectrum in hard X- 
rays by varying the coronal temperature, as shown in Fig llUI 
It is seen from this figure that the coronal temperature af- 
fects the hard X-ray spectrum, but does not change the soft 
X-ray continuum. In opposite, the change of Mj^t affects the 
form of only the soft X-ray spectrum (Fig llip and the hard 
X-ray is left unchanged, so we can infer the jet mass-loss 
rate Mj^t by adjusting the spectrum in the soft X-rays (the 
temperature of the jet base is assumed to be equal to that 
of the corona). Note that the change of both jet and corona 
parameters cannot compensate each other because the re- 
lated spectral components affect the resulting spectrum in 
different energy ranges. After that, we adjust the corona op- 
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Figure 11. The results of the simulations for different values of 
the jet mass loss rate: 3 ■ 10^^ g/s (the dotted line), 4 • 10^^ g/s 
(the solid line), 5 ■ 10^^ g/s (the dash-dotted line). The best fit 
is found for Mj^t = 4 ■ 10^^ g/s, Tcor = 0.24, Tcor = 19 keV, 
rcor = 6.4 ■ 10" cm. 
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Figure 12. The results of the simulations for different optical 
depth of the corona: 0.2 (the dotted line), 0.24 (the solid line), 
0.3 (the dash-dotted line). The best fit is found for Tcor = 0.24, 
Tcor = 19 keV, rcor = 6.4 ■ lO" cm, Mjct = 4 • lO^^ g/s. 



tical depth (FiE ll2p . The variations of coronal optical depth 
at Tcor < 1 do not change the form of the spectrum, but do 
change the coronal X-ray luminosity. So we fit the luminos- 
ity of the coronal component as shown in Fig |12l Finally, 
we constraint the accretion disk emission temperature at 
which the soft part of the model spectrum is still consistent 
with observations (Fig |13p . When the effective temperature 
of the accretion disk emission is too high, a soft X-ray excess 
appears, which is not present in observations. In Fig |14l we 
present the results of our simulations when the X-ray excess 
from the accretion disk partially compensates the low ener- 
getics of the jet. This figure clearly demonstrates that these 
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Figure 13. The results of the simulations for different values of 
the accretion disk effective temperature: 6.0 ■ 10^ K (the solid 
line), 9.4 ■ 10^ K (the dash-dotted line), 1.1 ■ 10^ K (the dotted 
line). The best fit is found for T^/ = 6.0 ■ lO'"^ K, Tcor = 0.24, 
Tcor = 19 keV, Tcor = 6.4 ■ 10" cm, Mj^t = 4 ■ lO" g/s. 
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Figure 14. Fitting the observational data for a model with low 
kinetic luminosity of the jet for different values of the accretion 
disk effective temperature: 9.410^ K (the dotted line), 1.1-10^ K 
(the dashed line). Tcor = 0.24, Tcor = 19 keV, rcor = 6.4-10" cm, 
Mjct = 3-1013 g/s. The solid line corresponds to Tcf = 6.0- 10^ K, 
Mjct = 4:- 10^^ g/s with all other parameters as shown for previous 
lines. 



model parameters do not give quantitative agreement with 
observations. 

The last three figures clearly show that the emerging X- 
ray continuum is very sensitive to slight changes of physical 
parameters and our results form an island in the parameter 
space. Further observations and more precise determination 
of the distance to the source, the component mass ratio, etc., 
can change the jet mass loss rate and corona density. 

We also should stress that the hard X-ray spectrum 
of SS433 by itself can be formally fitted in several physi- 



cal models, ranging from purely thermal emission from jets 
with changing temperature to variants of standard comp- 
tonization procedures from the XSPEC package. An exam- 
ple of such a descrip t ion w ith CompPS model is given in 
ICherepashchuk et al.l |20o3). However, when the physical 
parameters of the system are fixed by other observations, 
like those discussed in this paper, the description of hard X- 
ray spectrum by standard tools become not self-consistent. 
So the results of spectral modeling presented in this pa- 
per can be considered as a first attempt to find disk, jet 
and corona parameters in a self-consistent way. Clearly, fur- 
ther observations of the complicated variable emission from 
SS433 are needed to improve the model parameters. 

We would like to draw attention to some recent results 
in determining the binary compone nt masses in SS433. The 
optical spe ctroscopy of the system l|Hillwig and Gie^l2008l : 
iHillwig eF al. 2004) revealed the presence of absorption lines 
in the spectrum of the optical component identified as a 
~ ATI supergiant star. Observed orbital Doppler shifts of the 
absorption lines of the optical component allowed Hillwig & 
Gies to determine the mass ratio of the relativistic (Mx) 
and the optical (M„) components in SS433 q — M^/M^ ~ 
0.3 ± 0.11, implying the masses Mo; = 4.3 ± O.SAf© and 
Mv = 12.3 ± 3.3M0. These results are not the final ones 
(cf. the analysis by Cherepashchuk et al. (2008), who ob- 
tained Mx ~ 5M0), and the values of the binary compo- 
nent masses can change with different interpretation of op- 
tical spectroscop ic observations, as thos e authors mention 
in the conclusion iHillwig and Gied l|2008l ). Nevertheless, the 
smaller mass of the compact object compared to our choice 
{Mbh = 7M0) will result in decreasing of the correspond- 
ing Roche lobe size, and thus the accretion disk radius. But 
X-ray spectral properties of SS433 are determined by jets 
and corona above central parts of the accretion disk, so the 
X-ray spectrum of the source will hardly change significantly 
for smaller Mx- 

We conclude that the observed broadband X-ray con- 
tiimum of SS433 is dominated by thermal emission of the 
jet and rarefied hot corona around the inner parts of the 
supercritical accretion disk, with addition from the inverse 
Compton scattering on hot electrons of the corona. The jet 
mostly contributes in the soft X-ray band, while hard X-ray 
emission is dominated by thermal emission from the corona 
and inverse Compton scattering of soft thermal X-ray pho- 
tons. The emission from inner regions of the accretion disk 
also contributes to the spectrum, but mostly at energies 
below 1 keV. The observed broadband X-ray spectrum of 
SS433 is best fitted by the corona temperature ~ 19 keV 
with Thomson optical depth Tcor — 0.24 and the jet mass 
loss rate about 4 - 10^^ g/s which corresponds to the kinetic 
luminosity of the source 1.2 - 10''^ erg/s. 
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7 APPENDIX: MONTE-CARLO TECHNIQUES 

The Monte-Carlo method has another name: the statistical 
test method. An innumerable number of processes can be 
listed the outcome of which is not determined and has a 
probabilistic nature. When repeated, they will give us dif- 
ferent results. But these processes obey statistical laws that 
can be ascertained by repeating one process many times, i.e. 
by performing statistical tests. As the number of tests in- 
creases, the fluctuations decrease, and we can evaluate some 
value that we are interested in. For instance, the mean value 
of a random variable, its variance, the density probability, 
etc. The process of photon propagation through a medium 
has probabilistic character: the photon interacts with mat- 
ter, this interaction changes the photon properties (e.g. its 
frequency, the direction of motion) and one cannot predict 
exactly the outcome of the process, i.e. whether the pho- 
ton escapes from the medium, where the photon will es- 
cape, and its new frequency. Because of that the Monte- 
Carlo method is suitable for deriving the characteristics of 
the radiation field, such as its spatial and frequency dis- 
tribution. Of course, one can solve the radiation transfer 
equation to obtain these quantities using some approxima- 
tions for stationary and time- dependent problems, e.g. see 
ISunvaev and Titarchu^ l|l980l ), pTitar chuk (1988|), but in aU 
cases the geometry of the problem is rather simple and in 
many cases cannot be applied directly to sources with com- 
plec geometry of emitting/scattering regions, like SS433. 
Also, the finite-difference numerical solution of an integro- 
differential equation with frequency dependence is quite so- 
phisticated and requires much computational time. These 
are our reasons for choosing the Monte-Carlo method. 

One possible way to perform a single test with a photon 
is to simulate its trajectory in a medium. Initially a photon is 
placed somewhere within the computational domain, its fre- 
quency and the direction of motion are known, and we trace 
the photon as it propagates through the medium (below we 
shall discuss how one can do this). The moment when the 
photon leaves the computational domain, or its trajectory is 
terminated, is considered to be the end of the test and the 
longest time of one test is taken to be the duration of the 
Monte-Carlo experiment itself At. But there is a problem: 
any real source emits enormous amount of photons, so it will 
take correspondingly enormous amount of time to simulate 
all trajectories, and one will have to wait for years to get the 
result. It is clear that we should abandon real photons and 
introduce model ones. One possible way is to trace many 
photons combined in a group, i.e. a photon packet. At this 
point, another problem arises: when a pack of photons inter- 
acts with the medium, some photons are, for example, scat- 
tered, others are absorbed or produce an electron-positron 
pair in the field of an atomic nucleus, or something else. 
One can split the photon packet and follow the trajectory of 
each part, but there is another way out. In our simulations 
we used the indivisible photon packets method developed by 
ILucvI l|l999l ). According to this method, individual photons 
are grouped in the packets with constant energy 




Figure 15. Initial spectra: accretion disk emission (dotted line), 
jet free-free emission (dashed line), corona free-free emission 
(dash-dotted line). The solid line is the result of the Monte-Carlo 
simulation. Comptonization of emission on corona electrons is re- 
sponsible for 'hardening' of the spectrum. 



Such a trick helps us to avoid individual tracing of a large 
number of photons, especially in the low-energy part of the 
spectrum, since the number of photons in a packet increases 
with decreasing frequency. The process of interaction of a 
photon from the packet is chosen each time, according to the 
known probability of the variant. Such a method may look 
inexact, but the increasing of the number of photon packets 
makes less probable things happen, which approaches the 
result to the exact solution. 



7.1 The number of photon packets 

At the beginning of the Monte-Carlo experiment a number 
of the photon packets is placed within the computational 
domain. This number is arbitrary, but the higher this num- 
ber, the better accuracy can be achieved in the resulting 
spectrum. As a consequence of the statistical character of 
the experiment, the relative error in the result is inversely 
proportional to the square root of the photon packets num- 
ber, so to increase the accuracy by an order of magnitude 
one should increase the number of photon packets by two or- 
ders of magnitude. About 10* photon packets is enough to 
provide a smooth spectral curve that can be compared with 
observations. In our case there are photon packets of dif- 
ferent origin: free-free photons from the jet and corona and 
the accretion disk photons. Each 'photon type' has its own 
area of the initial spatial location and spectral distribution. 
In Fig llSl the initial spectra are presented together with the 
result of the Monte-Carlo simulation which is discussed in 
the main text. Each component makes its contribution to 
the source luminosity. The free-free emissivity of plasma in 
the jet or corona is 



eo ~ nhv. 



(10) 



dVnl{r) dn dukAT)B^{T) 



(11) 
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It should be taken into account that the temperature 
varies along the jet, so one should first integrate over the fre- 
quency and then over the jet volume. The corona is isother- 
mal and L ; f^^or can be calculated as a product of three inte- 
grals. The photon luminosity of the accretion disk was taken 
to be a fraction of the Eddington luminosity Ldisk = V^Edd 
with the factor 77 being inferred during the fitting procedure. 

One should set the number of photon packets that 
will participate in the Monte-Carlo experiment. After calcu- 
lating the total luminosity of all components 



Ltot — Ldisk + Lff,jet + Lff^cor, 

we determine the value eo / At as 



So _ Ltot 
At ^ lV~' 



(12) 



(13) 



where eo is the energy of the packet, At is the duration of 
the Monte-Carlo experiment. Once this value is calculated 
we can determine the number of photon packets representing 
diflterent spectral components as 

Li Li 



N, 



eo/At Ltot 



N. 



(14) 



The index 'i' stands for the 'type of emission' (corona or jet 
free- free or accretion disk emission). 



7.2 The distribution of photon packets over the 
frequency bins 

The frequency to the photon packet should be assigned. The 
necessary spectral range is divided into several frequency 
bins. Their number must be sufficiently high to provide a 
smooth spectrum, but small enough so that the average 
number of photons (packets) per bin Nph be larger than 
unity, Nph/Nun ^ 1. In our calculations we divided the 
spectrum into 500 equal bins on logarithmic scale. The ini- 
tial frequency was taken in the middle of the frequency bin. 
The number of photon packets in each frequency bin Nij is 
determined by the photon spectral distribution and is pro- 
portional to the contribution of the frequency bin to the 
emissivity of the spectral component: 



Nij/N, 



"3 00 



Li, I, dv 



(15) 



temperature distributions in the corona, the volume was di- 
vided into spherical layers with constant density and tem- 
perature. The corona was divided into 300 layers of constant 
thickness, and jet was divided by similar spherical cuts into 
about 10^ truncated cone layers. The contribution of each 
layer to the total emissivity was calculated according to 
(|11[) and a certain number of photon packets representing 
the free-free spectrum of each layer was introduced. In each 
layer its own free-free spectrum was maintained constant. 
The photon packets within a layer were distributed with 
equal probability for a free-free photon to appear in each 
point in the layer 

r = [(rf-r-tOli+^-f-i]^ 
cos6'i„ = (cos6ljet/2 — sin6'disfe/2) 72 + sin6'disfc/2 

ipin = 27r73 

Here and ri_i are outer radii of the i-th layer, 71,2,3 
are random numbers distributed uniformly over the interval 
(0, 1), 6 jet is the jet opening angle and Odisk is the accretion 
disk thickness angle. 

In the case of the accretion disk photon packets, the sit- 
uation does not differ drastically. They are emitted from the 
spherical segment with radius ro, and should be uniformly 
distributed over this surface: 




ro 

(1-s 
27r73 



9disfc 72)72 ~¥ sin 6 disk/ "2 



7.3.2 The initial direction of motion 

One should set the initial direction of motion for the photon 
packets. Free -free photon packets emitted in a plasma move 
isotropically (|Soborlll973t ). 

cos 9 — 274 — 1 
(fi — 27175 

In the Cartesian coordinates, unit vectors (orts) determining 
the initial direction of motion are written as 

( ex = sin 9 cos ip 
\ Sy = sin ^ sin ((3 



e, = cos f 



Initial photon packets from the accretion disk are assumed 
to move radially. 



All photon packets in each frequency bin are assigned the 
same frequency. During a Monte-Carlo experiment, the 
trial photons are chosen sequentially starting from the low- 
frequency bins to the high-frequency ones. 



7.3 Initial data for photon packets 

7. 3. 1 The spatial position 

Let us determine the spatial positions of the photon pack- 
ets. Obviously, in a homogeneous and isothermal medium 
the free-free photons should be distributed uniformly. But 
if the medium is inhomogeneous or/and non-isothermal, the 
situation is more complicated. One should divide the whole 
volume into regions of constant density and temperature. In 
the particular case of spherically symmetrical density and 



7.4 Simulation of the photon packet trajectory 

Simulation of propagation of photons through the medium is 
the essence of the Monte-Carlo experiment because at this 
stage the resulting angle and spectral distribution of the 
photon packets is formed. Once the initial conditions for a 
photon packet are set, one can find a randomly distributed 
optical depth the packet must travel before it interacts with 
matter using a simple Monte-Carlo formula 



T = — In 7, 



7 



is a random number. 



(16) 



Next, we should choose the velocity and the direction 
of motion of an electron that scatters the photon packet. 
For electrons with relativistic Maxwellian distribution, the 
velocity of electron can be calc ulated using the r andom 
number algorithm suggested by IPozdnvakov et al] l|l983h . 
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In the corona all directions of motion are equally probable, 
so the distribution of orts of the electron motion is isotropic 
(|7.3.2|l . The matter in the jet moves with the radial velocity 
0.26c. The distribution of directions of motion of electrons 
is isotropic in the local co-moving jet frame. So the resulting 
vector of motion of an electron is obtained by the Lorentz 
transformation of the isotropic ort (in the co-moving frame) 
to the rest frame. 

After that the coordinates of the photon-electron inter- 
action point can be found. We determine the total optical 
depth as 

r = y dln,{r)[(Tc + K{T)n,{r)] (17) 

For a given density distribution we can integrate this ex- 
pression analytically, taking into account that the photon 
packet trajectory is a straight line defined by the equation 

n+i = r i -f se„ (18) 

Thus we obtain a transcendental equation for the dis- 
tance s the photon travels between two subsequent interac- 
tions. This equation is solved numerically using the Newton 
method. After that we obtain the required coordinates of 
the interaction point from (|18p . We must check whether the 
interaction point lies within the region occupied by plasma 
(the computational domain). If it does not, then we trace 
the photon trajectory to the boundary of the computational 
domain, and collect all necessary information, such as pho- 
ton's frequency and escape angle. 

If the interaction point lies within the computational 
domain, we should choose the process of interaction of ra- 
diation with matter. We assume that the continuum X-ray 
spectrum of SS433 in the 3-90 keV range is formed due to 
the Compton scattering by free electrons and free-free ab- 
sorption/emission, so the choice should be done randomly 
between these two events. We select the physical event by 
means of a simple Monte-Carlo procedure. A random num- 
ber uniformly distributed over the interva l (0, 1) is selected, 
and if it satisfies the criterion (|Lucvll 19991 ) 



then the photon packet was considered to suffer a Compton 
scattering, otherwise it was considered to suffer a free-free 
absorption. If scattering occurs, we should determine the 
new f requency of the packet an d its new direction of motion 
as in IPozdnvakov et al] (|l977h . The energy of the photon 
packet should be changed with the change of frequency ac- 
cording to (|10p . so the number of photons in the packet 
remains constant. If the photon packet is absorbed, then it 
simply disappears, and we proceed to the next packet. 

This sequence should be repeated until the photon 
packet is absorbed or has left the computational domain. 
Then we gather all required information about the escaped 
photon packet. To complete the Monte-Carlo experiment, 
this procedure is applied to all photon packets. 
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